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FOREWORD 


This report describes developments in the analytical determination of unsteady supersonic 
aerodynamics for a class of interfering wings using a triangular element representation of wings and 
diaphragms. ..... 
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DEVELOPMENT AND APPLICATIONS OF 
SUPERSONIC UNSTEADY CONSISTENT 
AERODYNAMICS FOR INTERFERING 
PARALLEL WINGS 

by 

Kari Appa * and G. C. C. Smith + 

Bell Aerospace Company 
Division of Textron Inc. 


SUMMARY 


The analytical development of unsteady supersonic aerodynamic influence coefficients, 
(AIC’s), for isolated and nearly parallel interfering coplanar and noncoplanar wings is described. 
Numerical formulations based on triangular discretizations of wings and diaphragms are handled in 
a kinematically consistent manner. Examples of isolated wing cases are compared with respect to 
AIC’s and flutter boundaries. AIC’s for interfering wings are compared where corresponding results 
are available. 

Computer programs for the interfering case are described in a companion User’s Manual, 
NASA CR-112184, and presented in a Programmer’s Manual, NASA CR-1 12185. 


INTRODUCTION 


The evaluation of aeroelastic effects requires accurate determination of oscillatory aerody- 
namic, stiffness and inertia properties of a flight vehicle. Stiffness and inertia properties of geometri- 
cally complex structures may be computed to a high degree of accuracy using modern finite element 
methods. The determination of unsteady aerodynamic forces in various flow regimes on the other 
hand, involves more complex differential equations, and their integral solution equivalents. Satisfac- 
tory solutions, analytical for simple planforms, and numerical for more complex planforms have 
evolved for isolated planar wings. Much effort has been devoted to more complex situations involv- 
ing intersecting surfaces, interfering surfaces, dihedralled surfaces, and wing-body interference. The 
summary papers of References 1 and 2 attest to the intensive interest and efforts applied. 

Garrick and Rubinow (Ref 3) pioneered the development of solutions to such boundary value 
problems by the superposition of elementary solutions of the nature of pulsating acoustic sources 
distributed on the boundaries, from which have developed many techniques discussed in the literature 
References 1 and 2 indicate that integration of the differential/integral equations for practical and in- 
creasingly complicated wing and wing-body configurations is a significant problem even with present 
day computational facilities. Therefore, numerical integration methods are continually being modi- 
fied and improved. 


* Chief, Dynamic Analysis 


+ Chief Engineer, Structural Dynamics 


In supersonic flow, numerical integration procedures have been based on the representation of 
domains of dependence in terms of discrete “elements” of various types in each of which the down- 
wash has generally been approximated as a constant. 

These types include 

- square boxes (e.g., Pines, Dugundji and Neuringer, (Ref 4)) in which the sides are parallel 
and normal to the free stream 

- characteristic boxes (e.g., Smith (Ref 5) and Zartarian, Hsu and Voss (Ref 6)) in which 
parallelogram elements have their sides parallel to the Mach lines 

- rectangular (“Mach”) boxes (e.g., Li, Ref 7) such that the box diagonals are parallel to 
the Mach lines. 

Zartarian and Hsu (Ref 8) made extensive studies on the application of rectangular, character- 
istic and Mach box methods, and discussed various working rules for the application of the Mach box 
method. 

Ashley (Ref 9) described the application of superposition methods to interfering surfaces and 
the use of Mach box methods for intersecting planes. 

Moore and Andrew (Ref 1 0) and Donato and Huhn (Ref 1 1 ) developed computer programs to 
calculate the velocity potentials and generalized forces for isolated, interfering and intersecting planes. 
Stark (Ref 12) contributed very significantly to the evaluation of downwash terms along a subsonic 
leading edge in the characteristic box approach. 

The general approaches listed above have at least two obvious factors leading to inaccuracies 
in numerical computations 

- the planform and associated diaphragm boundaries cannot be well represented with 
elements tied to characteristic or flow directions 

- the approximation of constant downwash within an element. 

While both inaccuracies can be alleviated by the use of a larger number of smaller elements, 
it is obvious that this is at the expense of a large increase in the number of equations to be handled, 
and will result in rapidly increasing computational effort. 

Appa and Smith in References 13 and 14 have discussed a different approach to planform and 
diaphragm discretization and numerical procedures which have the following advantages over the 
methods referred to above: 

- constancy of the wing grid system at all Mach numbers 

- continuous distribution of downwash terms across element boundaries 

- well-represented wing and diaphragm boundaries 

- flexibility in the choice of diaphragm elements 
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- consistent aerodynamic forces via the principle of virtual work 

- convenient aeroelastic formulation for flutter synthesis problems (e.g., optimization) 

The following sections relate to analytical developments and computer applications for such 
a method in the case of isolated and interfering nearly parallel coplanar and noncoplanar wings at 
supersonic Mach numbers. 

Results for unsteady aerodynamic forces from the present method are compared with pub- 
lished data from other methods. 

Additionally, the method is applied to flutter of planar isolated and interfering wings and 
compared (for isolated wings) to experimental results. 

The general intent of the work is to illustrate the potential for more economic and accurate 
determination of supersonic unsteady aerodynamics particularly for complex configurations. 

An operational computer program which was developed in the course of the work has been 
used on IBM 360/65 and CDC6600. 
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SYMBOLS 


A', B\ A, B 
a 


a u> a C 
b u , bg 


c 

F 

I 



K 

e 

M 

M, 

n 

n 

N 

N 

P 

P 


Element and system pressure matrices 

A Boolean matrix describing topological connections of the elements in the 
structure 

Speed of sound in free stream 

Slopes of the upper and lower sides of a triangle 

Constants in the equations for upper and lower sides of a triangle 

Semichord (as in Ref. 21) 

Wing chord, Figures 3, 4, 5 
Kernel, defined by equation (5) 

Integrals; also unit matrix 

Reduced frequency 

kM 2 

P 2 

Stiffness matrix 

Reference length 

Mass matrix 

Mach number 

Unit normal vector 

Number of degrees of freedom 

Number of terms in interpolation or modal series 

Summation variable 

Pressure 

Nodal force vector 
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SYMBOLS (CONT) 


q 

q 

o 

Q, Q 
s 

r, s 
t 

T 

U 

w 

W 

w 


x, y,z 
x, z 

X = (x, y, z) 
Z (x, y, t) 

Z (x, y) 
a 


P 

r, r 
8 

V 

X 


Displacement vector 

Dynamic pressure, q* a reference pressure = Vi p* a* 2 M 2 
Generalized aerodynamic coefficients, discrete and interpolated, respectively 
Surface of integration 
Hyperbolic radii, Eqn (6) 

Dimensional time 
Transformation matrix 
Flight speed 
Downwash velocity 

Induced downwash influence coefficient matrix 
Virtual work 

Dimensionless rectangular coordinates 
x and z coordinates of apex of tailplane (see Fig. 9) 

Position vector of receiving point 

Displacement in space-time, nondimensional, referred to £ 

Displacement amplitude (nondimensional) 

Damping rate, (see X) 

i / 2 

= (M 0 2 -l) 

Transformation Matrices Eqns (19), (38) 

= q 0 */q 0 , dynamic pressure ratio 

A column vector of generalized coordinates 
= a + icu 


X 


XC a£ co£ 

~ = ~ + i“ , a complex reduced frequency 
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SYMBOLS (CONT) 


M 


M. v 


s = rj, n 


p 

/O 

a 

a 

•P 

0 

# 

X 

* 

* 

CO 


CO 


Ratio of structural mass to air mass 
Defined by Eqn. (A7) 

Dimensionless running coordinates in the fore Mach cone 
Position vector of source in the fore Mach cone 
Density of air 

Kinematic displacement vector of an element 
Source strength (scalar) 

A column vector of source strength at nodes of elements, and the complete system 
Velocity Potential - a scalar 
Column Vector of Velocity potential 

Element V.I.C. Matrix relating velocity potential to displacement for isolated wings, 
and to source strength for interfering wings (see Eqns. 1 1 and A-5 respectively) 

Transformation Matrix of Natural Modes 

Defined in Eqn. (A-9) 

System source strength to system displacement transformation matrix 
Circular frequency, radians per second 

as defined in Ref. 22 
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n 

Subscripts 

cr 

i. j 
£,m 
n, r 
Re, Im 


An interpolation matrix 
Critical 

Denote elements or summation variables 
Denote modal numbers 
Summation variables 
Real and imaginary parts 
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SYMBOLS (CONT) 


R 

Reference value 

s 

Wake 

2, u 

Lower and upper limits 

V 

Derivative with respect to r? 

Superscripts 


d 

Diaphragm 

L,U 

Lower, Upper 

t 

Transpose 

w 

Wing 

♦ 

Sea level value (indent on p, q Q , a 0 ) 

Notation 

i) 

A column vector 

t ) 

A matrix 

r j 

A diagonal matrix 

[ 1* 

Transpose of a matrix 

[ r 

Inverse of a matrix 

A( ) 

Difference between upper and lower value 

Note: 

In general terms, a lower case symbol represents a continuous variable, a bold-faced 
symbol with subscript/superscript represents an element nodal value matrix, a 
bold-faced symbol without subscript represents a “system” matrix (i.e., that result- 
ing from assembly of all element nodal matrices). This usage applies to p, q, w, 
p, a, <p, <//, in particular. 
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DEVELOPMENT OF THE METHOD 


Detailed derivation and some applications of “finite element” idealizations to the determina- 
tion of unsteady aerodynamic coefficients appear in References 13 and 14. Application to isolated 
and interfering planar wings is given in this report. The basic elements used to represent the planform 
of a wing in a discretized form may be quadrilateral and triangular elements. However, the simpler 
triangular element forms the basis in this work. The elements can be of different sizes and orientations 
to “best fit” both wings (independent of Mach number) and diaphragms. 

In finite element analysis, functionals within an element are expressed in terms of nodal (dis- 
crete) values and derivatives at the nodes by means of interpolation functions. The accuracy of the 
analysis depends on the degree to which the properties of continuity of the functionals across the 
element boundaries are satisfied. This is related to the order of appearance of the functional and its 
derivatives in the integral equation. Since no derivatives of source functionals are involved in the 
determination of potential, source strength distribution is taken to be expressed linearly within the 
elements in terms of their nodal values. 

Since modal data are also frequently expressed as discrete nodal displacements, the virtual 
work integrals determined in this work use linearly interpolated potential and modal data to give 
linearly consistent generalized forces. 


Velocity Potential 


The linearized equations of motion in unsteady supersonic flow referenced to a rectangular 
coordinate system, Figure 1 , are given by 


n m a-> , 32 ^ _ ^2 8 V . _J_ iV 

( ' ° } 3x 2 3y 2 3z 2 a Q 3x3t a Q 2 3t 2 


( 1 ) 


where <p, M Q and a Q represent the perturbation velocity potential, Mach number and speed of sound 
respectively. Solutions to equation (1) are sought for harmonic motion of the lifting surface. In 
many common flutter solution procedures the eigenvalue solutions subsequently determined using 
assumed real reduced frequency aerodynamics are complex (implying diverging and decaying oscilla- 
tions) for which the aerodynamics are inconsistent. This suggests, (compare Richardson, Reference 15) 
the formal determination of aerodynamics for complex exponential motions, which, though not carried 
through to flutter solutions herein, requires motion of the lifting surface in the form 

Z (x, y, t) = Z (x, y)e Xt 


where X = (a + ioo), (3) 

and a is a damping factor, and u> the circular frequency of oscillation. (See also Ref. 1 6). 

The deflection Z is related to the source strength through boundary conditions discussed later 
(see for example, Eqn. 13). 
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The solution for the velocity potential following Reference 3 is given by 

V (X, t) = - ^ JfF (X, 2) o (2) d £ dr? e Xt 
S 

where F (X, 2) = 



( 4 ) 


(5) 


is the kernel function, 

r = (x - £) - 0 [ (y - t}) 2 + (z-f) 1 ]^ 

s = (x-$) + 0 [(y-T?) 2 + (z - f )* ] 1/2 

are hyperbolic radii and a (2) is the nondimensional source strength distribution. The symbols X = 

(x, y, z) and 2 = (£, r), f) denote respectively the vector positions (dimensionless with respect to 
£) at which the velocity potential is evaluated and the influence source strength is situated. The sur- 
face integration in equation (4) covers the lifting surface and associated diaphragm regions in the fore 
Mach cone emerging from X. Generally, closed form integration of (4) is not possible, and numerical 
integration techniques using various “box” schemes are employed. In the current analysis, the integra- 
tion area is divided into a number of triangular elements (Figure 1). Summing over all elements the 
velocity potential within an element i may be written 

^,0 = ~~ J // F(X i 2 j ) ■ a (2j) • d £ d 77 • e Xt (7) 

j Sj 

where the subscripts i and j denote receiving element i and influencing element j respectively. By 
means of interpolation functions fi (2) the source strength distribution a (2j) within an influencing 
element can be expressed as a linear combination of its nodal values (Tj as 

a (2j) = S2 (2) 0 j (8) 

where S2 , (see Eqn. A-2) is a row matrix of interpolation functions and oj is a column vector 
of nodal values. The nodal values of the source strength for each element can be abstracted from a 
system source vector cr by 

a j = aj a (9) 

in which aj is a Boolean Matrix (Ref 1 7) related to the topological assembly of the elements. Sum- 
ming over all the elements in the fore Mach cone region, the velocity potentials at the nodes of the 
receiving elements are given by 

0i =-Ufi X a j c (10) 

j 



where a is a nondimensional nodal source vector and 

=ir// F(X n<S)dS 

S; 


( 11 ) 
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Figure 1. Layout of Grid System 







is the velocity potential influence coefficient matrix for the ith element nodes due to the source dis- 
tribution in the jth element. Numerical integration of $ jj is described in the Appendix. Then con- 
sidering all the elements, the velocity potential for the system can be written as 

0 = -Ufi ($1 { <7 } (12) 

where 0 is a column vector of nodal velocity potentials and $ is the total velocity potential influence 
coefficient matrix. 

Boundary Conditions 

The source strength nodal vector <7 in Eq. (12) must be such that the following boundary 
conditions are satisfied : 

• the kinematic boundary condition on the lifting surface given by 

dip DZ 

1? - * dT (13) 

where 

dp 

—— is the derivative of the velocity potential normal to the surface at the field point X, 
on 

and 

_D_ = _9_ U _3_ 

Dt " 3t + £ 3x ° 4) 

is the material derivative in the streamwise direction. 

• in the case of a subsonic leading edge, the velocity potential difference in the diaphragm 

A <t > d = Supper - *lower>a must be zero 05) 

• the pressure difference in the wake 

A P S = Supper ' Plower>s must be zero. (16) 

The boundary conditions for interfering wings require special considerations described later. 
Subsonic Leading Edge Case 


Antisymmetry of the nodal source distribution cr and of 0 with respect to the plane of the 
wing in the lifting case implies that Equations (12) and (13) may be Used to relate potential and source 
strength jumps in partitioned matrix notation as 


1 A K j = . T jo [fww *wd 1 j A ^w 

' A ^d f L $ dw ^dd J i A<7 d 
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where the subscripts ‘w’ and ‘d’ denote the nodes situated on the wing and diaphragm respectively. 
Since in the diaphragm A 0 d = 0 (Ref. 18), the source strength vector A o d from the second of Equation 
(17) is 

{A<r d } = -r{Aa w } (18) 

where 


r = 1 $ 

dd dw 


(19) 


Using this solution in the first set of matrix Equation (1 7), the velocity potential distribution on the 
lifting surface is given by 

A0 W =-^3^ A<r w (20) 

in which 

^ww" ^ww ' ^wd ^ (21) 


Kinematic Boundary Conditions 
The Isolated Wing Case 

__ It is shown in Ref. 3 that the source strength distribution, a e^ is proportional to the downwash 
DZ At 

2 — = we . Let/*; denote a column vector of nodal kinematic deformation modes of an element 
Dt ' J 

j in the fore Mach cone region. Then the interpolated deformation within the element j can be ex- 
pressed in terms of the element nodal values as 


Z($j,Tjj,t) = ft (£, V)fj e 


At 


( 22 ) 


If the mean position of the lifting surface is parallel to the x-y plane, the downwash from Eq. (13) 
is 


1 dtp v 0 U 8 — 

2 0? 0t 2 0£ 

or in nondimensional form using Equation (22) 

— J - (Aft + ft^)/ 5 j 
- a2 ia>2 

in which A =(— + -jj— ) is complex reduced frequency, and 
0ft 

ftt . 

« 0 ? 

Thus for the isolated case, the source strength distribution ae^* in Equation (4) will be replaced by 


(23) 


(24) 


(25) 


v (Hj) = (A ft + ftj) Pj 


(26) 
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Utilizing the relation between element nodal displacements/^ and system displacements q as 
given by the Boolean matrix aj 

/ 5 j =a j q (27) 

the system velocity potential analogous to Eq. (12) may be expressed in matrix form as 

A<j> = -2Ufi [<fr] q (28) 


Interfering Nearly Parallel Wings 


The determination of source strengths in the case of interfering planes is more complex than 
in the isolated case. Ashley in Reference 9 has discussed in detail some practical assumptions for 
relating the source strength distribution to the kinematic boundary conditions. In the present analysis 
the following assumptions based on Reference 9 are made: 

(i) convenient diaphragms can be established between the lifting surfaces and the Mach 
hyperbolae 


(ii) induced normal wash at one diaphragm due to the presence of the second lifting surface 
can be disregarded. 


Considering two planar interfering wings (Figure 2 ), the source distribution oj w on the lower sur- 
U 

face of wing 1 and a 2 W on the upper surface of wing 2 are unaffected by the interference. Then 
these source strengths using Eq. 23 can be expressed as 


L £ DZ, 

° lw " 'u dT 


(29) 


U __fi_ DZj 
° 2w ~ U Dt 


(30) 


The source strength distributions on the upper surface of wing 1 and the lower surface of wing 2, be- 
cause of interference, are modified by the induced downwash in their respective planes. The induced 
downwash on the upper surface of wing 1 due to source distributions on the lower surface of wing 2 
and its diaphragm may be written as 

w 12 = W 12w g 2w + W 12d a 2d ( 31 ) 

where 

Wj 2 is a nodal column vector of induced downwash on wing 1 and its diaphragm 

Wj 2 is the velocity influence coefficient matrix for the wing 1 due to interference of source 
distributions on the lower surface of wing 2 and its diaphragm. 
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Figure 2. 


i-ayout of Grid System - 2 Wings 
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The subscripts w and d denote control points on the wing and its associated diaphragms, and 
the superscripts U and L denote upper and lower surfaces, respectively. Derivation and numerical 
integration of the downwash influence coefficient matrix Wjj for the triangular element is described 
in the Appendix (see Eqs. A13-A20). 

In the following work, in accord with assumption (ii), the effect of induced downwash wj 2 
in the diaphragm has been disregarded. Since the wing planes are assumed to be nearly parallel, the 
side wash contribution to the computation of w j 2 has also been neglected. Then the source strength 
required to satisfy the total normal velocity condition on the upper surface of wing 1 is given by 

U £ DZ, L L 

a lw " u "dT ' W 12w°2w- W 12d *26 ( 32 > 


Similarly the required source strength distribution on the lower surface of wing 2 is given by 


L r fi DZ 2 U U I 

2w |_U Dt 21w CTlw " 21d °ld J 

Eqs. (32) and (33) can be rearranged as 


- “ 

1 u 

/ £ DZ! 




U 

(ll W 12w 

a lw 

I U Dt 


0 - W 12d 


°ld 




+ 





/ L 

j^_ dz 2 




L 

; W 21w l 22 

( °2w 

' U Dt 


_ W 21d 0 


°2d 


(33) 


(34) 


The relation between the source strengths in the diaphragm and wing is given by Equation (18). Re- 
writing this equation as 


U 

L 




u 

L I 

a ld 

- °ld 


r 


a lw ‘ 

a lw / 

U 

L 




U 

L ( 

a 2d 

' a 2d , 


- - 


a 2w ' 

°2w ) 


and using the condition of downwash continuity in the diaphragm, i.e., 

U L 

°d = ’ °d (36) 


the required source strength in the diaphragm may be written as 


U 

a ld 


- -] 

/ u 
CT lw 


- - 


L 

°lw 

L 

a 2d 


r 

/ L 

( °2w 


f 


U 

°2w 


(37) 
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where 


1 r-i °"i r i 1 0 

= - r 

2 Lo i J L J Lo -i . 


( 38 ) 


Eqs. (29) and (30) then give 


\ 

u 1 

- - 


U ) 

- 

°ld ( 



°lw / 


) = 

f 


\ + 

r 

L \ 





^2d ) 

- 


°2w ) 

- - 


_£_ H± 

U Dt 

_£ DZ 2| 

U Dt 

» » l 

U L 

Substituting for ctjj and a 2 $ > n Eq. (34), the following matrix equation is obtained: 


-W 


w 


12w 


21 w ! 22 


-W 


1 2d 


W 21d 0 


U 

°lw 

L 

°2w 


£ DZ, 


U Dt 

-£ DZ 2 
"Tj Dt~ 


r- n 




-o 

<N 

sT 

i 

O 


r 


L W 21d 0 . 


- - 



£ DZ, 


U Dt 

-£ Hi 

U Dt 


£ DZ, 


(39) 


(40) 


X U£-,\ 

The input downwash terms are calculated for the system structural displacements 

U Dt 

3q 


q at each grid point, and the slopes — , by 

ox 


£ Hk — -v f l . j 

u'Bf ‘ k|q(+ IT 


(41) 


Using Eq. (41) in Eqs. (29), (30) and (40), the difference in the source distribution across the wing 
surfaces can be written as 

Aa w J = [ *] { q \ (42) 

where 'i' is a transformation matrix relating source strength and displacement vectors. The velocity 
potential differences across the wing surfaces are then obtained by substitution of Eq. (42) into Eq. 
( 20 ). 


Concerning Subsonic Trailing Edges 

In this case the source strength distribution in the wake can be determined from the pressure 
continuity condition across the wake sheet and the trailing edge potentials. It may then be expressed 
in terms of the distribution on the wing and diaphragm and condensation of velocity potential in the 
wake performed similar to Equation (21). 
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For interfering cases, induced downwash matrices can similarly be modified. However, the 
subsonic trailing edge case has not been included in the present analyses arid programs. 

Consistent Nodal Forces 

The velocity potential within the i^ element can be expressed in terms of nodal values 
A$j and the interpolation matrix S2 as 

Ay> (Xj) = [ O 1 { } e Xt (43) 

Since the pressure distribution within the i 1 * 1 element from the linearized Bernoulli’s equation 
is given by 

f 3 U 3 T 

_] (44, 

from Equations (43), (23), (24), the pressure distribution is 

P (Xj) = + -y- [In + n x ] A 0 j (45) 

The virtual work done on the i** 1 element by the pressure through virtual displacement fiSZ (Xj) 
using Eq. (22) is 

5Wjj = pU^-Spj 1 //[!«* ft + fttft*] dxdy A^ (46) 

S i 

Summing over all elements on the lifting surface and using Eq. (28), the virtual work principle yields 
the nodal forces, 


p = pUB[AA + B] A <(> 

.= -2pU 2 £ 2 [1 A + B ] ( ] q 


where A and B, real pressure matrices, are assemblies of element matrices defined by 

A - S a n A n a n 
n 

B X a n B n a n 
n 


The element matrices are given by 




/ /ft* ft ds 
s j 


(47) 

(48) 

(49) 

(50) 
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//«*«* ds 


and ‘n’ represents summation over wing elements only. 


Generalized Forces 


If qg and q m denote two nondimensional nodal deformation modes, the virtual work done 
in a virtual displacement 5qg is, using Eq. (47) 


5W Cm - C6qg p r 


= -2p U 2 C 3 5qg [ X A + B ] [ $ ] q r 
- -2 p U 2 C 3 Sq* Qg m q m say. 


The matrix of dimensionless generalized force coefficients is then defined by 

1 d / ~ 

2p U 2 C 3 3q fi [ 3q m J 


the elements of Q. 
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FLUTTER ANALYSIS 


The flutter equations of motion for a flight structure may be written 
Kq + Mq + 4 (^fi A(k,M 0 )q = 0 

where K, M, A are (nxn) stiffness, inertia and aerodynamic matrices respectively, 
q is a (nxl) displacement vector 


(54) 


£ is a reference length 

2 3 

pU 2 pa o Mo 

qo=V = — 

k =-^y- is the reduced frequency 


is the dynamic pressure and 


p is the air density 

U is the true flight speed. 

For small oscillations, the deformation of the structure can be approximated by a linear 
combination of a few of the lower frequency natural modes as 


q = X?7 


(55) 


where X is a nxN modal matrix, (N«n) and rj is a vector of generalized coordinates. Taking the 
motion to be of the form 

_ Xt 

ri = T?e (56) 

Eqn (54) becomes, for zero damping at the flutter boundary (a = 0), 

K iiVi - w 2 + 4q 0 £ 2 Qy (k, M^j = 0 (57) 

for i = T. . . N, where Ky, My and Qy are the generalized stiffness, inertia and aerodynamic matrix 
coefficients respectively. 


Using a reference frequency ccr and a reference dynamic pressure q* = — — | — ^2. t 

where * parameters refer to sea level condition, say, Eqn (57) can be rewritten in nondimensional 
form as 





[Q(k, M 0 )J 



V 


(58) 


in which 8 = -^2. 

<13 


is a dynamic pressure ratio, the ccj are the natural frequencies of the K - M system 



and co, v are the eigenvalue and eigenvector of the system respectively. Since the evaluation of Q 
requires a priori knowledge of the flutter frequency, the eigenvalue problem cannot be solved 
directly. Several iteration schemes have appeared in the literature, the method for the determination 
of flutter eigensolutions employed in this work being as follows: 


For a given Mach number, compute Q(k, M 0 ) for a range of reduced frequencies, k. Since 
the elements Qjj vary smoothly with k, each Qy can be approximated by a power series/ spline fit 
over a chosen range of k, as 


Qij (k) 




(59) 


For a given Mach Number, the variables in Eqn. (58) are the density p and the speed of sound a Q . 
Generally, the relation between p and a 0 is given in the form of standard atmospheric tables or 
tunnel operating characteristics, effectively reducing the unknowns to one variable, say p. Rewrite 
the definitions of dynamic pressure and reduced frequency in terms of Mach number as 


and 


q Q ="JPU 2 =-jpM 0 2 a 0 2 '= pga*»M Q 2 g 

2 

, _ co£ _ <o£ 

U M 0 a 0 


(60) 


(61) 


The iterative flutter solution in the atmospheric flight condition begins with the assumption 
of a low density p and the corresponding speed of sound a^. Assuming some starting value of co = 
ojr, k and 5 are then determinable. For this k, Qy’s are interpolated. With these quantities all 
the eigenvalues of Eqn. (58) are determined. Since Q is complex, the eigenvalues are generally 
complex, with their imaginary parts positive for decaying motion, which is inconsistent with the 
harmonic motion assumed. In order to determine consistent results, density is now increased and 
the above process repeated until one or more of the imaginary parts of co change sign from positive 
to negative. The reduced frequency ‘k’ is then recomputed using the smallest real part of co which 
has a negative imaginary part. The density p is then interpolated between those values of p which 
yield positive and negative imaginary parts of to to obtain acceptable convergence to zero of the 
imaginary part of ‘to’. In practice rapid convergence to a real reduced frequency is obtained since 
the flutter formulation given by Eqn. (58) is generally not sensitive to small changes in k. In most 
cases linearly interpolated values of p are adequate for determination of the critical value. The 
critical dynamic pressure q ocr , critical flight speed U cr , and the flutter frequency cof are finally 
determined. A typical flow chart is presented in Figure 13. 


Presentation of the flutter boundary in the current examples is in terms of the conventional 
V 1 

stiffness-altitude parameter, -*— = where p, the ratio of structural mass to air mass is chosen by 

the user. ^ 
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RESULTS AND DISCUSSION 


The numerical results presented in this section have been obtained from two independent pro- 
grams developed for isolated and interfering cases. 

The characteristics of and parameters related to the various results are summarized in Table 1 . 

Isolated Wings 


Aerodynamics for Isolated Wings 

The accuracy of any unsteady aerodynamic method is critically illustrated by comparison of 
velocity potentials. Figures 3 to 5 compare the velocity potential distribution for a cropped delta 
wing with a low subsonic-leading edge at M 0 = 1 .054 for steady and oscillatory cases. Five consistent 
triangular elements compare very well with 30-box results from Reference 1 2. For the same cropped 
delta wing, generalized aerodynamic coefficients are presented for heave and pitch about a mid-chord 
axis at Mach numbers 1.118 and 1 .201, in Table 2 and Figure 6. For such low Mach numbers and so 
few elements, the generalized coefficients deduced from the present method are in encouragingly good 
agreement with Reference 12 with the exception of moments due to pitching. 

Table 3 shows the generalized coefficients for a rectangular wing of aspect ratio 2.0 at M 0 = 

2.0 in heave and pitch motion with k = 0.3, and in steady flow at M 0 = 1 .2. These have been compared 
with References 19 and 1 1 , respectively. Vector plots of generalized forces are presented in Figure 7 
and show generally good agreement. 

The next example considers an AGARD swept wing of aspect ratio 1.45. The wing motions 
considered are heave, pitching about center chord, quadratic bending in the chordwise direction and 
quadratic bending in the spanwise direction, at two Mach numbers and two frequencies in each case. 
Comparisons have been made with Reference 20 which uses 30 to 40 characteristic boxes and Refer- 
ence 19 which uses 17 chordwise Mach boxes whereas the present method uses only six triangles. 

Table 4 shows the generalized coefficients for M 0 = 1 .2, k = 0.5. For the first three modes good agree- 
ment has been obtained with References 20 and 19 while the fourth mode compares well with Refer- 
ence 20 only. It appears that the Mach box results in Reference 1 9 are wrong or inaccurate for this 
mode. Table 5 shows Qy : s at M 0 = 1.2, k = 1.0 for the same modes. Once again good agreement with 
Reference 20 is obtained. 

Tables 6 and 7 show generalized coefficients at M 0 = 2.0 for k = 0.5 and k = 1 .4 respectively. 
Comparisons have been made again with References 19 and 20 in the former and with Reference 19 
in the latter case. Results in Table 6 are in good agreement with Reference 20 but for the fourth mode 
(spanwise bending), again Reference 19 agrees with neither. 

Flutter Solutions for Isolated Wings (NASA “HT-7” Example) 

The basic purpose in the development of consistent aerodynamic forces is ultimately to im- 
prove the accuracy of aeroelastic solutions. Generalized forces have been used to determine the flut- 
ter boundary for a wing of aspect ratio 2.50 and taper ratio 0.3, known at NASA Langley as the HT-7 
configuration. The leading edge sweepback angle is 50.5°. Mode shapes, generalized masses and 
natural frequencies are given in Reference 21. Flutter boundaries at M Q = 1.3, 1.573 and 1.64 were 
determined using the flutter equation (58). The altitude stiffness ratio b 0 co 2 J~jT and the flutter 
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TABLE 1 

SUMMARY OF EXAMPLES 
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TABLE 2 

COMPARISON OF AERODYNAMIC GENERALIZED COEFFICIENTS Qy 
FOR A CROPPED DELTA WING 
M 0 = 1.118, 1.201, k = 1.05 (3 =0.3) 

Modes: q, = 2/7, q 2 = x-x midchord 


Method 

Ref. 12 

Present 

No. . of <L 





Elements 

30 

6 

M 0 

a»i 

Real 

Imaginary 

Real 

Imaginary 


On 

0.0602 

0,3163 

0.0519 

0.2849 


Qi j 

1.1603 

-0.0915 

1.0811 

-0.0452 

1.118 






Q 21 

-0.0504 

-0.0746 

-0.0353 

-0.0390 


q 22 

-0.1631 

0.3923 

-0.0578 

0.3039 


Qu 

0.0646 

0.2992 

0.0560 

0.2763 


Q, 2 

1 .0727 

-0.1159 

1.0134 

-0.0678 

1.201 






Q 2 , 

-0.0437 

-0.0738 

-0.0359 

-0.0497 


q 22 

-0.1574 

0.3684 

-0.0847 

0.3132 




TABLE 3 

COMPARISON OF AERODYNAMIC GENERALIZED COEFFICIENTS Qy 
FOR RECTANGULAR WING OF ASPECT RATIO 2 
M 0 = 2.0, k = 0.3 and M 0 = 1.2, k = 0 
Modes: qj = 1.0, q 2 = (x - c/2) 
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TABLE 4 

COMPARISON OF AERODYNAMIC GENERALIZED COEFFICIENTS 
Qjj FOR AGARD SWEPT WING OF ASPECT RATIO 1 .45 
M 0 = 1.2, k = 0.5 

Modes: q! = 1.0, q,= (x-c r /2), q 3 = (x-c r /2) 2 , q 4 = y 2 


Method 


No. of £ 
Elements 


Present 

Mach Box Ref. 19 

— » 

6 

17 

30+ 

Real 

Imaginary 

Real 

Imaginary 

Real 

Imaginary 

-0.0862 

3.10 

0.0119 

3.473 

-0.0228 

3.3506 

3.46 

2.18 

3.801 

1.681 

3.671 

1.7325 

3.6 

-1.69 

3.532 

-1.185 

3.353 

-1.9538 

0.0357 

0.766 

-0.0894 

0.702 

0.0319 

0.7478 

-0.32 

0.13 

-0.2695 

-0.016 

-0.2832 

0.0083 

0.37 

2.53 

0.227 

2.475 

0.244 

2.459 

2.58 

-1.11 

2.852 

-1.055 

2.703 

-1.0507 

-0.0359 

0.108 

-0.0983 

0.022 

-0.0216 

0.0764 

-0.194 

0.5525 

-0.1474 

0.589 

-0.164 

0.5532 

0.68 

1.45 

0.7242 

1.276 . 

0.6824 

1.305 

1.4 

-0.312 

1.444 

-0.226 

1.357 

-0.2489 

0.0272 

0.0985 

-0.0692 

0.064 

-0.0182 

0.0827 

-0.1135 

0.723 

-0.3301 

-0.388 

-0.0223 

0.8380 

0.82 

1.00 

\ 

-0.3376 

1.627 

0.9066 

0.574 

1.1 

-0.386 

0.6071 

-0.151 

0.9925 

-0.1059 

-0.0432 

0.232 

-0.0978 

-0.111 

-0.018 

0.2528 





















TABLE 5 

COMPARISON OF AERODYNAMIC GENERALIZED COEFFICIENTS 
Qjj FOR AGARD SWEPT WING OF ASPECT RATIO 1.45 
M 0 = 1.2, k = 1.0 

Modes: q, = 1.0, q 3 = (x-c r /2), q 3 = (x-c r /2) 2 , q 4 = y 2 


Method 

Present 

Char. Box Ref. 20 

No. of <£ 





Elements 

6 

30+ 

Q ii 

Real 

Imaginary 

Real 

Imaginary 

Qn 

-0.565 

3.52 

-0.4042 

3.5477 

Ql2 

4.12 

1.804 

4.0734 

1.5832 

Q.3 

2.6 

-0.081 

2.56717 

0.15633 

Ql4 

-0.159 

0.785 

-0.13087 

0.72916 

Q 2 i 

-0.724 

0.735 

-0.70852 

0.50164 

Qj2 

0.988 

1.684 . 

0.8165 

1.73847 

Q 23 

2.082 

0.192 

2.21076 

0.09960 

Ql4 

-0.216 

0.212 

-0.17584 

0.13608 

Q 3 1 

-0.376 

0.888 

-0.33974 

0.80848 

Q 3 I 

0.962 

0.93 

0.91061 

0.8632 

Q3 3 

1.238 

0.454 

1.19235 

0.43602 

Q 34 

-0.153 

0.18 

-0.12188 

0.12959 

CLl 

-0.322 

0.960 

-0.145 

0.89538 

Q 42 

0.719 

0.808 

0.97682 

0.52965 

Q 43 

0.808 

0.139 

0.82489 

0.23027 

Q 44 

-0.196 

0.297 

-0.13382 

0.26639 




TABLE 6 

COMPARISON OF AERODYNAMIC GENERALIZED COEFFICIENTS Qj- 
FOR AGARD SWEPT WING OF ASPECT RATIO 1 .45 
M o = 2.0,k = 0.5, Modes: q, = 1.0, q 2 = (x - c r /2), q 3 = (x - c r /2) 2 , q 4 = y 2 


Method 

Present 

Mach Box 

Ref. 19 

Characteristic Box Ref. 18 

No. of <t 







Elements 

6 


17+ 

40 

°ii 

Real 

Imaginary 

Real 

Imaginary 

Real 

Imaginary 

Qn 

0.0698 

2.6047 

0.06007 

2.565 

0.05999 

2.52563 

i Qj 2 

2.6094 

0.4950 

2.572 

0.5374 

2.53299 

0.51324 

0,3 

1.5535 

0.6932 

1.566 

0.6912 

1.51857 

0.65485 

Q 14 

-0.0121 

0.6462 

-0.01934 

0.5779 

-0.01374 

0.57410 

Q 2 i 

-0.0057 

0.4767 

-0.001025 

0.4388 

-0.00307 

0.44491 

^22 

0.4913 

0.7677 

0.45666 

0.7695 

0.46181 

0.75181 

Q 2 3 

1.4937 

0.2990 

1.533 

0.2829 

1.48154 

0.27395 

Q 24 

-0.0230 

0.1655 

-0.02136 

0.01302 

-0.01853 

0.14258 

Q 3 i 

-0.0095 

0.5410 

-0.00299 

0.5293 

-0.00676 

0.50594 

°32 

0.5482 

0.4688 

0.5422 

0.4457 

0.51692 

0.44825 

Q 3 3 

0.8621 

0.3366 

0.8699 

0.3102 

0.84474 

0.30134 

Q 3 4 

-0.0234 

0.1133 

-0.01891 

0.1028 

-0.01730 

0.10143 

Q41 

0.0186 

0.7421 

-0.1180 

-0.1415 

0.01719 

0.67412 

Q 42 

0.7319 

0.2103 

-0.1348 

0.6712 

0.66558 

0.19623 

Q4 3 

0.5752 

0.3510 

0.3856 

0.07844 

0.53698 

0.32057 

Q 44 

-0.0224 

0.2793 

-0.05113 

-0.03105 

-0.01956 

0.23862 
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TABLE 7 

COMPARISON OF AERODYNAMIC GENERALIZED COEFFICIENTS 
Qj FOR AGARD SWEPT WING OF ASPECT RATIO 1 .45 
M 0 = 2.0, k = 1.4 

Modes: q t = 1.0, q, = (x-c r /2), q, = (x-c r /2) 2 , q A = y 2 


Method 

Present 

Mach Box Ref. 19 

No. of £ 





Elements 

6 

17 

°ii 

Real 

Imaginary 

Real 

Imaginary 

Qn 

-0.060 

2.4 

-0.1529 

2.35 

Qj 3 

2.550 

0.77 

2.5010 

0.8112 

Q.3 

1.511 

0.556 

1.480 

0.5726 

0.4 

-0.134 

0.666 

-0.1862 

0.6014 

Oil 

-0.336 

0.484 

-0.3306 

0.4262 

Q 33 

0.625 

0.856 

0.5708 

0.8512 

Qj3 

1.356 

0.28 

1.357 

0.2858 

Q 34 

-0.179 

0.218 

-0.1774 

0.1791 

Q 3 . 

-0.258 

0.574 

-0.2473 

0.541 1 

Q 33 

0.647 

0.517 

0.6291 

0.4906 

Q 3 3 

0.766 

0.334 

0.7360 

0.3294 

Q 34 

-0.167 

0.172 

-0.1496 

0.1413 

04. 

-0.488 

0.682 

-0.7363 

0.1362 

04, 

0.665 

0.328 

0.1201 

0.5232 

0,3 

0.625 

0.271 

0.2283 

0.1796 

Q 44 

-0.155 

0.332 

-0.3423 

0.00649 













Nondimensional Velocity Potential <j>/ bU Q • Real Part 



Figure 4. Real Part of Velocity Potential on a Cropped Delta Wing In Translation 
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Nondimensional Velocity Potential <j)/b\J Q - Imaginary Part 







Imaginary (Qj|3 Imaginary (Qjj) 

M 1.0 r 



(a) IV^ =2.0, k = 0.3 (b) M Q = 1.2, k = 0.0 


Figure. 7- Vector Plot of Aerodynamic Generalized Coefficients Qy for 
a Rectangular Wing of Aspect Ratio 2 
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to reference frequency ratio ( w/w 2 ) are shown in Figure 8. Solutions obtained from supersonic 
kernel function and piston theories (Reference 21) are also shown for subsonic leading edge and 
supersonic leading edge cases respectively. The flutter boundary from the kernel function theory at 
a Mach number of 1 .573 yields a conservative result while the piston theory is unconservative. The 
flutter boundary obtained from the consistent approach lies in between these two results, encourag- 
ingly close to the experimental flutter point. In the subsonic leading edge case, the present method 
gives a frequency ratio very close to the kernel function approach. In the supersonic leading edge case 
it suggests a transition toward piston theory solution. The experimental flutter frequency is lower 
than from the present solution. 


Interfering Wings 

Two sets of computations were made using the program for interfering wings. The first 
example was that considered by Woodcock and York in Reference 22. (See Figure 9.) The second 
(NASA example) was a combination of two tapered wings of aspect ratio 6.8 and leading edge sweep 
of 16.0°. (see Figure 10.) No experimental comparisons were available for this latter case. However, 
flutter calculations have been carried out to illustrate the effect of separation distance at two Mach 
numbers. Additionally, the effect of static deformation on flutter has been roughly assessed at one 
Mach number. 

Generalized Aerodynamics (Interfering Wing-Tail Example) 

The reference 22 wing-tail combination is shown in Figure 9. The aspect ratios of wing and 
tail are 2.31 and 4.62 respectively, the tail span being half the wing span. The apex of the tail is 
situated at x = 2 N /lTand z = 0.5, which corresponds to d = 0.866 and z = 0.5 in Reference 22. The 
number of elements used in Reference 22 and in the present analysis are respectively, 300 on the wing 
with 50 on the tail, and 49 on the wing with 16 on the tail. 

Computations were made at a Mach number of 1 .44 for which the wing leading edge is sub- 
sonic and the tail leading edge is supersonic. The modes are heave and pitch of the whole configura- 
tion and heave and pitch of the tail alone, the pitch axis in both cases being located at 2/3 of wing 
root chord. The computation was performed at a reduced frequency of k = 0.01 as in Reference 22. 

In order to give a rough indication of some of the interference effects, a few heave damping 
coefficients for the wing alone and tail alone are quoted from the present work, and References 22 and 
23. 


Results are presented in Table 8 together with corresponding figures from Reference 22, and 
some comparable results from Reference 24 forz = 0. 

No objective comments are offered at this stage except that generally comparable results are 
obtained for coefficients of significance. 

Flutter of Interfering NASA Wings 

The two wings considered are of identical planform with zero stagger (See Figure 10). Aspect 
and taper ratios are 6.8 and 0.364 respectively, and the leading edge sweep is 16°. 

Flutter calculations were performed at Mach numbers of 1 .45 and 1 .60 for various separation 
distances including no interference. The results are presented in two forms on Figures 1 1 and 12, as 
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1.0 1.2 1.4 1.6 1.8 2.0 

Mach Number 


Figure 8. Flutter Boundary in Terms of the Stiffness - Altitude Parameter and Ratio of Flutter 
Frequency to the Second Natural Frequency versus M for HT-7 Wing 
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Figure 10. NASA 1A and 2A Configuration and M 0 =1.6 Aero Grid 


















flutter speed against separation for two Mach numbers, and flutter speed against Mach number for 
various separation distances. Additionally at M 0 = 1 .60 only, the effect on flutter speeds of deflec- 
tion due to thickness loading was assessed roughly. This was done by first determining the flutter 
dynamic pressure for zero thickness. At this dynamic pressure the steady interference loading due to 
thickness was determined and applied to the wings to determine their approximate deflected equilib- 
rium positions. Flutter calculations were then repeated for these deflected positions by recomput- 
ing zero-thickness generalized forces, etc. 

Modal data for unit generalized masses of the two wings are given in Tables 9A and 9B for the 
first four modes of each wing, these being the modes used in the flutter calculation. 

Figure 1 1 shows the effect of separation on reduced flutter speed at M 0 = 1 .45 and 1 .60, for 
which asymptotic flutter speeds without interference are 1 .77 and 1 .87 respectively. At M 0 = 1 .45, 
interference commences at the root at a dimensionless separation of 0.95 and has spread to the tip at 
a separation of 0.348. At M 0 = 1 .60, the separations are 0.80 for the root and 0.293 for the tip 
interference. 

The results for M 0 = 1 .60 suggest that the effect of interference, which starts at the root, is 
first small, as might be expected. As the inteference effect reaches beyond mid-span, the flutter speed 
falls much more rapidly, the rate of fall decreasing again as the interference effect spreads out to the 
tip. The same general conclusions hold for M 0 = 1 .45. However, at this Mach number the total reduc- 
tion in flutter speed is more drastic and is probably affected by the shortcomings of linearized theory 
which predicts negative two-dimensional aerodynamic damping for a range of reduced frequency be- 
low M 0 =/T The wings in the present example are not too far removed from the unswept two-dim- 
ensional case. 

Figure 1 2 crossplots flutter speed against Mach number and also sketches the asymptotic non- 
interference boundary. 

Both Figures 1 1 and 12 include the effect on flutter speed at M 0 = 1 .6 of the static interfer- 
ence. The effect is small and the single iteration used for static deflection is probably adequate. This 
effect has not been determined for M 0 = 1 .45. 

An interesting feature of the flutter modal vector behavior was noticed as the separation 
decreased. Without interference, the thinner, lower frequency Wing 1 A fluttered at a lower 
reduced speed than Wing 2A (i.e., at 1.87 and 1.77 for Mach numbers of 1.60 and 1.45 respectively 
-see Figure 1 1). The flutter mode was dominated by Mode 3, its first torsion mode. 

Z 

However, for the smallest separation case ( — = 0.214) of the Mach 1.60 calculations, the 

2b 

flutter mode and frequency were dominated by the thicker wing 2A. The flutter frequency was 
2837 rps compared to the zero airspeed Wing 2A Mode 3 (torsion) frequency of 2819 rps. 

The significant contributions to the modal vector were: 

Wing 1 A - Mode 3 (first torsion) +0.080 

Wing 2A - Mode 1 (first bending) - 0.083 

Mode 3 (first torsion) + 0.928. 


TABLE 9A 

MODE SHAPES AND FREQUENCIES FOR NASA WING MODEL 1 A 


Mode 

1 

2 

3 

4 

Frequency 

346.0 .rps 

1290.0. rps 

2121.4 rps 

3169.0 rps 

Grid Point 

- 

- 

- 

- 

1 

0.0 

0.0 

0.0 

0.0 

2 

-0.248374E 00 

0.846053E 00 

0.340652 E 00 

0.156274E 01 

3 

0.0 

0.0 

0.0 

0.0 

4 

-0.661 762 E 00 

-0. 132569E 01 

-0.1 62827 E 02 

-0.955371 E 01 

5 

-0.148843E 01 

-0.574492E 01 

-0.501 730E 02 

-0.306777 E 02 

6 

0.373210E 01 

-0.1 10920E 02 

0.256854E 02 

-0.238349E 02 

7 

0.5 101 67 E 01 

-0.1 51 07 2 E 02 

0.166748E 02 

-0.290629E 02 

8 

0.600098 E 01 

-0.188923 E 02 

-0.3831 27 E 01 

-0.338759E 02 

9 

0.714894E 01 

-0.260340E 02 

-0.358998E 02 

-0.482695E 02 

10 

0.751287E 01 

-0.329570E 02 

-0.7 52936 E 02 

-0.650456E 02 

11 

0.144421E 02 

-0.37711 IE 02 

0.664367 E 02 

-0.624758E 02 

12 

0.168970E 02 

-0.433420E 02 

0.358953E 02 

-0.6251 67 E 02 

13 

0.189125E 02 

-0.494378E 02 

-0.280577 E 01 

-0.635001 E 02 

14 

0.202884E 02 

-0.557970E 02 

-0.480677 E 02 

-0.659036E 02 

15 

0.210618E 02 

-0.629261 E 02 

-0.969772E 02 

-0.717563E 02 

16 

0.327125E 02 

-0.6566 15E 02 

0.1 07401 E 03 

-0.631332E 02 

17 

0.356679E 02 

•0.688935 E 02 

0.59481 IE 02 

-0.503386E 02 

18 

0.3801 85E 02 

-0.728230E 02 

0.537 178E 01 

-0.3791 78E 02 

19 

0.397843E 02 

-0.779427 E 02 

-0.5308 14E 02 

-0.274941 E 02 

20 

0.409142E 02 

-0.844757 E 02 

-0.1 13283E 03 

-0.197513E 02 

21 

0.5841 88E 02 

-0.713794E 02 

0.1421 17E 03 

-0.106175E 01 

22 

0.614552E 02 

-0.701302E 02 

0.804984 E 02 

0.196408E 02 

23 

0.639868E 02 

-0.705903E 02 

0.145253E 02 

0.4024 10E 02 

24 

0.659998E 02 

-0.730238E 02 

-0.541 276E 02 

0.602483E 02 

25 

0.674978E 02 

-0.777540E 02 

-0.123549E 03 

0.794476E 02 

26 

0.91 2201 E 02 

-0.31 0397 E 02 

0.158395E 03 

0.598651 E 02 

27 

0. 940856 E 02 

-0.257317E 02 

0.883828E 02 

0.758822E 02 

28 

0.965729E 02 

-0.226300E 02 

0.157768E 02 

0.941429E 02 

29 

0.986789E 02 

-0.219322E 02 

-0.581 623E 02 

0.1 14858E 03 

30 

0.100418E 03 

-0.237869E 02 

-0.132125E 03 

0.1 38293 E 03 

31 

0.129880E 03 

0.676545E 02 

0. 144149E 03 

-0.820722E 01 

32 

0.1 32357 E 03 

0.748086E 02 

0.730007 E 02 

-0.304004E 01 

33 

0.134593E 03 

0. 7 987 29 E 02 

0.785800E 00 

0.718163E 01 

34 

0.136588E 03 

0.827491 E 02 

-0.717218E 02 

0.229209E 02 

35 

0.138360E 03 

0.834467 E 02 

-0.143789E 03 

0.443335 E 02 

36 

0.171953E 03 

0.210059E 03 

0.989585E 02 

-0.242098E 03 

37 

0.173924E 03 

0.21 6496 E 03 

0.345747 E 02 

-0.241725E 03 

38 

0.175761E 03 

0.221376E 03 

-0.298829E 02 

-0. 235332 E 03 

39 

0.177465E 03 

0.224694E 03 

-0.941 044E 02 

-0.222813E 03 

40 

0.179052E 03 

0.226521 E 03 

-0.157824E 03 

-0.204166E 03 




TABLE 9B 

MODE SHAPES AND FREQUENCIES FOR NASA WING MODEL 2A 


Mode 

1 

2 

3 

4 

Frequency 

476.6 rps 

1716.0 rps 

2819 rps 

4211 rps 

Grid Points 

- 

- 

- 

- 

1 

0.0 

0.0 

0.0 

0.0 

2 

-0.216321 E 00 

0.733357 E 00 

0.293900E 00 

0.1 3436 IE 01 

3 

0.0 

0.0 

0.0 

0.0 

4 

-0.5764 53E 00 

-0.1 15600E 01 

-0.141411E 02 

-0.832933 E 01 

5 

-0.129031 E 01 

-0.500474E 01 

-0.435684E 02 

-0.266567 E 02 

6 

0.323887 E 01 

-0.962833E 01 

0.22281 5E 02 

-0.2071 37 E 02 

7 

0.442801 E 01 

-0.131 154E 02 

0.1 44487 E 02 

-0.252496E 02 

8 

0.520807 E 01 

-0. 163985E 02 

-0.333866E 01 

-0.294 104E 02 

9 

0.620506E 01 

-0.225997 E 02 

-0.31 161 IE 02 

-0.418901 E 02 

10 

0.652201 E 01 

-0.2861 14E 02 

•0.653320E 02 

-0.564257 E 02 

11 

0.125306E 02 

-0.327 159E 02 

0.5761 16E 02 

-0.542299E 02 

12 

0.146600E 02 

-0.375993 E 02 

0.31 1273E 02 

-0.542434 E 02 

13 

0.164083E 02 

-0.428873E 02 

-0.244496E 01 

-0.550728E 02 

14 

0.176038E 02 

-0.4841 08E 02 

-0.417061 E 02 

-0.571416E 02 

15 

0.182764E 02 

-0.5460 15E 02 

-0.841242E 02 

-0.621955E 02 

16 

0.283770E 02 

-0.569429E 02 

0.931 467 E 02 

-0.547628E 02 

17 

0.309410E 02 

-0.597486E 02 

0.515929E 02 

-0.436474E 02 

18 

0.3298 10E 02 

-0.631616E 02 

0.4663 18E 01 

-0.328564 E 02 

19 

0.3451 37 E 02 

-0.676064E 02 

-0.4603 18E 02 

-0.237896E 02 

20 

0.354956E 02 

-0.732784E 02 

-0.982352E 02 

-0.170476E 02 

21 

0.5067 14E 02 

-0.6 1887 IE 02 

0.1 23267 E 03 

-0.918818E 00 

22 

0.533058E 02 

-0.608072 E 02 

0.698350E 02 

0.170479E 02 

23 

0. 555023 E 02 

-0.612102E 02 

0.126208E 02 

0.3493 13E 02 

24 

0.572493E 02 

-0.633248 E 02 

-0.469 159E 02 

0.523054E 02 

25 

0.585499E 02 

-0.674306 E 02 

-0.1071 10E 03 

0.689820E 02 

26 

0.791 174E 02 

•0.268896E 02 

0.137385E 03 

0.518804E 02 

27 

0.816033E 02 

-0.222904E 02 

0.766725E 02 

0.657801 E 02 

28 

0.837612E 02 

-0.1 9605 IE 02 

0.137075E 02 

0.81 6333 E 02 

29 

0.855884E 02 

-0.1 90043 E 02 

-0.5041 15E 02 

0.996200E 02 

30 

0.870982 E 02 

-0.206159E 02 

-0.1 14543E 03 

0.1 19965E 03 

31 

0.1 12643E 03 

0.587043E 02 

0.1 25004 E 03 

-0.721 169E 01 

32 

0.1 14791E 03 

0.649023 E 02 

0.63301 2E 02 

-0.270577E 01 

33 

0.1 16731E 03 

0.692881 E 02 

0.672285E 00 

0.618898E 01 

34 

0.1 18462E 03 

0.71 7767 E 02 

-0.662088E 02 

0.1 98692 E 02 

35 

0.11 9999 E 03 

0.723772E 02 

-0.1 24703 E 03 

0.384641 E 02 

36 

0.149126E 03 

0.182181E 03 

0.857724E 02 

-0.21001 IE 03 

37 

0.1 50836 E 03 

0.187757E 03 

0.299216E 02 

-0.209643E 03 

38 

0.152429E 03 

0.191982E 03 

-0. 259948 E 02 

-0.204055E 03 

39 

0.1 53908 E 03 

0.194855E 03 

-0.817059E 02 

-0.193157E 03 

40 

0.1 55284 E 03 

0.1 96436 E 03 

-0.136979E 03 

-0.1 76952 E 03 



Reduced Flutter Speed 



Figure 1 1. Variation of Reduced Flutter Speed with Vertical 
Separation at M Q = 1.45 and M Q =1.6 
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CONCLUSIONS AND RECOMMENDATIONS 


Further applications of an improved consistent approach to the derivation of unsteady 
aerodynamic forces for planar wings have been illustrated, and an operational computer program 
produced. Similar techniques have been developed for nearly parallel interfering planar wings, and 
an operational program has been used for the study of flutter characteristics of a pair of wings as 
affected by separation and steady airload deflection effects. 

In general terms, it appears that the approximations made in satisfying boundary conditions 
for interfering wings are warranted for the limited range of applications considered, but further 
evidence of their validity should be obtained as more comparative numerical examples become avail- 
able. 


The encouraging comparisons with References 11,12, 19 and 20 confirm general conclusions 
concerning the advantages of the present approach in reducing the numerical size of matrix equations 
necessary for acceptable accuracy. 

With respect to the excellent experimental flutter correlation obtained in Figure 8, the ex- 
ample stands alone and requires further applications before general conclusions are possible. The 
same remarks apply concerning the Kernel Function - Piston Theory transition noted in the text. 

The interference algorithm has illustrated the probability of strong detrimental effects of 
proximity on the flutter of interfering wings. 

The flutter algorithm used offers potential for more direct computer operator interaction 
than more common approaches and is worthy of additional development. 

The present programs are not optimally designed from the computer time point of view, but 
have generally evolved with the developments reported. No computer time comparisons with 
Reference 22 are available, for example, but again, the marked difference in grid sizes is apparent. 

The expansion of the present approach to include more complex interacting and interfering 
configurations, subsonic trailing edges and interfering wakes is feasible. It is probable that the 
present approach will be even more advantageous from the numerical point of view in these more 
complex circumstances. The basic method could additionally be reformulated in terms of an inter- 
grated potential approach (References 25-27) avoiding the use of diaphragms. 

An iterative form of static aeroelastic solution for deforming, interacting surfaces is feasible. 
In these problems, the domain of dependence of a particular point may alter significantly with the 
“static” deformation, and require an iterative type solution to the resulting equation. Flutter 
solutions could continue to be determined as linearized problems about the static equilibrium 
position. 


APPENDIX 


Derivation and numerical integration of the velocity potential and downwash matrices for 
basic triangular elements are described in this section. Details of the programming for the solution 
of the total problem are contained in Reference 28. 


Velocity Potential - Isolated Wing 


For thin wings in supersonic flow, the nondimensional source distribution due to perturbed 
motion from Eq. (26) is given by 


o (S j) = 


w (5 ,) 
U 


L'xft(S) +n 5 ] p j 


(A- 1) 


where the interpolation matrix ft (a j) for the triangular element is 


n (S) = [ $ rj > 1 ] T 




*1 

Vi 1 

in which 

T = 

£2 

ha 1 



£ 3 

173 1 


3 ft 

r 


and ft ^ = 

3 ? 

•= [ 

1 0 0 


T 


(A-2) 


(A-3) 
( A-4) 


Using the source displacement relation given by Eqn. (A- 1 ) in Eqn. (4), the velocity potential in- 
fluence coefficient matrix relating nodal source strengths to nodal displacements can be written 
similar to Eqn. (1 1) as 




£ x ft + ft 



dS 


( A-5) 


(Note the distinctive definitions of for isolated and interfering wings). 


Integration of Eqn. (A-5) is performed numerically using Gaussian 5-point quadrature 
coefficients. The elements cut by Mach lines emerging from the position Xj were truncated by a 
small width and integrated analytically using characteristic coordinates. More detailed information 
is contained in Reference 13. 


Velocity Potential - Interfering Wings 


For interfering wings the source distribution is not known until the kinematic boundary 
conditions as given by Eqn. (13) are satisfied. This requires the normal derivative of the velocity 
potential. When this is performed the kernel F in Eqn. (5) will exhibit a higher order singularity 
which usually requires more elaborate numerical integration. However, this situation can be avoided 
by using the equality credited to Watson (References 8 and 1 1), 

A 


cos 




a_ 

dv 



( A-6) 


\/M 2 - v 2 
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where m 2 = (x - £) 2 - 0 2 (z - D 2 


(A-7) 


v = & (y-r?) 


kM 2 

A K1Vt O 
k = j3 2 


k x i / v \ 

and p (M, v) = J c (“ M ) sin ( 

+ J 2r ( sin f 2rsin ‘ 1 (f)) 


r - 1 r 


in which J 2 r is a Bessel function of even order. 

Using (A-6) in Eq. (5) the velocity potential influence coefficient matrix from Eqn. (11) 
is given by 

1 r u a V r u dip a „ (A- 10) 

4>ij =— f e - ik < x i-£j) / )"5T 3,?d ^ 

T?g 


Substituting from (A7) and (A8) and integrating by parts, the velocity potential becomes 

*u A 

4>ij = ~ J c lk ,x rSj»[s2 (|,r| u )lKtt! u )-na,Dli) * (Ml) (A-ll) 

% % 


“'tin j 


The upper and lower limits of r? are given by 

^u = a u ^ + b u (A- 12) 

n = a C ? + bg 

where a u and ag are slopes of the upper and lower boundary lines, (i.e. line 1-3 and line 1-2 res- 
pectively) and b u and bg are the corresponding constants (See Figure 14 A). The initial position of 
each element is assumed to be the average of its z coordinates if the wing is not wholly parallel to 
the x - y plane. The integrand in Eq. (A 1 1) is analytic in the fore Mach cone region and Gaussian 
type numerical integrations have been employed in all cases. 
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C. Partial Element Integration 


Figure 14. Integration Scheme for Triangular Elements 





Downwash Integrals 


Assuming that the lifting surface is approximately parallel to the x-y plane, the normal 
velocity at the surface is given by 


\t 1 dip e 

w e = T IT 


Xt 


(A- 13) 


Using Watson’s relation (Eqn. (A6))in Eqn. (4), and performing the differentiation with respect to 
z, the nondimensional downwash matrix can be Written as 

« - 


L 


V u 

Vo 


W(t,V) dr? d£ 
3 z 


(A- 14) 


The derivatives in (A 14) can be changed to different bases more convenient for integration 
purposes. From (A7) 


3i p = 

b\p 

JiL 


3z 

bH 

3z 

(A 15) 

b\p_ 

b\p 

bn 



bn 


(A 16) 


and therefore 

P 2 (z-n J±_ 

dz (x-$) 3£ (A 17) 

Substituting Eqn. (A 15- 17) in Eqn. (A 14), the downwash integral, after performing integration by 
parts and simplifying, becomes 

= -3< z -» f I, - I, 1 ( A18 > 

7T 

£u A 

f „ -ik (x - £) ( 

JS — (^a,T? C ) -^tt,TJ u )) (A- 19) 

(x-S) 

+ ( n tt, T7g ) t ( i, - n (£, T? u ) * a, t? u )M 

(X - 5) ’ 


w ij 


where 1 t = 
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and I 2 = SI 


^ui 

/ 

n 


e ~^ (x ~*u> (£ u , I?) - e " ik(x '¥._ * 


«£, V) ' 


(* - $ u ) 


(X - £g) 


/ 


* (£, h) e lk (X ~ * } /(l +ik(x-£)\ d£ 


(x-£) 


(x-£) 


dr? 


(A-20) 


The first two terms in I, vanish when taken on the boundary of the triangle. The limits of integ- 
ration in I 2 are as shown in Figure 14A, B, and C. Since the singularities have been eliminated 
in these integrals the limits of integration for the elements can extend up to the Mach cone. 

Numerical integration for Eqns (A 11), (1 9) and (20) in the £ and tj senses proceed towards the element 
side furthest from the respective axis in all cases (See Figures 14A, B). Gaussian integration is per- 
formed for the now-analytical formulation of elements totally within and cut by the Mach cones. 

The limits of integration of the latter are determined by the conditions r or s = 0 as defined in 
Eq. (6). These limits are indicated on Figure 14C. 
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